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Abstract: This paper deals with the numerical analysis of two one-way systems 
derived from the general complete modeling proposed by M.V. De Hoop. The 
main goal of this work is to compare two different formulations in which a 
correcting term allows to improve the amplitude of the numerical solution. It 
comes out that even if the two systems are equivalent from a theoretical point 
of view, nothing of the kind is as far as the numerical simulation is concerned. 
Herein a numerical analysis is performed to show that as long as the propagation 
medium is smooth, both the models are equivalent but it is no more the case 
when the medium is associated to a quite strongly discontinuous velocity. 

Key-words: Microlocal numerical analysis, one-way formulation, acoustic 
waves 
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Propagation one-way a amplitude prservee dans 
des milieux heterogenes 

Resume : Dans ce travail, on s'intsse a deux formulations one-way de I'equation 
des ondes acoustique qui ont cte construites a partir du modele complet one- 
way propose par M.V. Dc Hoop. Lc principal objectif de cette etude est de 
comparer les deux formulations dans IcsqucUes on a introduit un terme permet- 
tant d'ameliorer le calcul de I'amplitude de la solution numerique. II ressort 
de I'analyse que meme si les deux systemes sont equivalents du point de vue 
theorique, il n'en est rien au niveau des performances numeriques. On montre 
en particulier que tant que le milieu de propagation est regulier, les deux modeles 
se comportent identiquemcnt et que les differences sont nettes si le milieu com- 
porte des heterogeites. On pent done en conclure que le precision de I'amplitude 
est tres sensible a la formulation du modele. 

Mots-cles : Analyse micro-locale numerique, formulation one-way, ondes 
acoustiques 
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1 Introduction 

The numerical solution of seismic acoustic waves propagation in heterogeneous 
media is generally based on the solution of the second-order full wave equation. 
Then the second-order wave equation can be solved completely by using a fi- 
nite difference scheme and it is well-known that such an approach results in a 
high computational burden, especially in the case of three-dimensional prob- 
lems. This surely explains that a lot of people prefer to solve an approximate 
problem which involves either a truncate expansion of the solution [lU| or an 
approximate wave equation arising from the factoring of the exact one [HI [13] . 
In the simplest case of a homogeneous medium, the solution can be obtained 
from the inversion of a system of one-way wave equations [I] . It is a very inter- 
esting way of solving the wave equation because it is based on the decomposition 
of the wave into a down-going part and an up-going one which reproduces the 
physical phenomenon very faithfully. During several years, people (see [7j and 
its references) tried to extend this approach by introducing correcting terms in 
the model to account for heterogeneities into the propagation medium. In 1996, 
M.V. De Hoop [5] derived a new formulation based on the micro-local analysis 
which allowed to derive a complete system of one-way wave equations coupling 
by exact correcting terms. The use of the theory of pseudo-differential opera- 
tors makes the derivation of the system easy. Nevertheless, its plain writing is 
complicated because it involves the composition of pseudo-differential operators 
which means that each term is defined from an asymptotic expansion. Thus 
in practice the complete system is approximated by truncating the asymptotic 
expansions. Such an approach may look as if it complicates the solution as 
compared to the now well-controlled solution of the full wave equation. But its 
formulation allows to unpack the multiples from the primary refiections which 
is of outstanding importance for the geological interpretations and allows to 
reduce the computation time. Paper [S] has been followed by numerous publi- 
cations and among them, we refer to as |16j in which one can find a complete 
bibliography on the topic. As far as the numerical solution is concerned, it is 
associated to the inversion of an approximate system generally based on only 
keeping the main term in each asymptotic development. J. Le Rousseau was 
the first to obtain accurate snapshots (see [TS] and [1]). However if using his 
numerical method for the computation of arrival times, one gets erroneous re- 
sults on the amplitudes level. More recently, Zhang et al. [HI have proposed 
a corrected one-way wave equation which allows to compute the correct ampli- 
tude of the acoustic pressure. This new one-way wave equation is obtained from 
the factorization of the full wave problem. The equation is factorized by using 
a WKBJ solution in which the amplitude is taken into account as well as its 
phase. Their new formulation of the one-way system includes a correction term. 
In this work, we intend to show that the amplitude of the numerical solution 
can be corrected by adding a transmission term in the system, proposed in [15j . 
The heterogeneities of the medium are modeled as discontinuities of the veloc- 
ity which is supposed to vary in all the directions. The correcting term can be 
included in the system by two ways and we show that one of them is optimal. 
In fact we show that the best improvement of the amplitude is obtained when 
including the transmission operator into the right-hand-side of the system. This 
might come upon the reader since the model which is the nearest of the one of 
Zhang and al. [18^ is obtained by including the transmission operator into the 
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one-way equations, i.e. into the left-hand-side of the system. But some numeri- 
cal tests indicate that the two approaches are close in the case of smooth media. 
The paper is divided into 7 sections, plus this one and a concluding part. The 
next one deals with the initial model whose unknown is the acoustic pressure 
and its transformation into a first-order system by introducing the velocity as an 
unknown also. The third part is devoted to the reduced system whose derivation 
is based on selecting the depth variable as the leading direction and by plug- 
ging the other terms into the frequencies domain after using a Fourier-Laplace 
transform. The fourth part concerns the first-order approximation of the re- 
duced system. By accounting for the complete coupling terms of order 0, we get 
two equivalent systems which are interesting to consider since their numerical 
solution can be obtained by two different ways. The fourth next parts deal with 
some numerical aspects which are essential for the method. We have chosen to 
neglect the description of the propagation because we intend to focus on the 
transmission operator. The numerical tests are developed in the 2D case but 
we mention that some 3D test have been performed in |16| . 

In the following, we use standard notations for the micro-local analysis of 
classical pseudo-differential operators and we refer to [l^ for their definitions. 
We only precise that the symbol of an operator P is denoted by a{P) and its 
principal symbol is ap{P). 

2 Initial model 

The analysis of the waves propagation is an efficient tool for imaging the soil. 
Assume the region of interest H. is located between the surface of the earth 
{z — 0} and a given depth {z = Zmax}- The phenomenon of propagation 
is supposed to be generated by a source located at the top of fi. Then the 
discontinuities of the propagation velocity can be defined by computing the 
time arrivals which are recorded by a set of receivers located at a given depth. 
By assuming that Vl is surrounded by two homogeneous regions Vlsup and rJm/ 
(see Fig[T]), the receivers do not record any wave propagating below or above f2. 
Hence only J7 is under study. 

The phenomenon is governed by the wave equation set in IB?: 



-^dfp - div{Vp) = S dans lR^x]0,Tl 
p(0, x) = dtpiO, x) = dans M^, 



(1) 



whose solution is the acoustic pressure p — p{t,x). The nonnegative variable 
t denotes the time while x —* {x,y,z) stands for the position vector. In the 
following, x' —* (x, y) designates the transverse variable. The propagation 
velocity v :— v{x) is constant-valued, equal to c*"^ in flsup and Cinf in flmf. 
In fi, V := c{x) varies in all the directions. The source S* is a Ricker function 
modelling an impulsion, 

9^ ( --^^(t-uf 



S{t,x)^S{x-Xs)—(e-'' 



where 5x^ denotes the Dirac distribution at Xs —^ {xs,ys,0), a — ttz^* and 
i, = l/z^*. The constant z/* is a given frequency. We use standard notations for 



INRIA 



True amplitude one-way propagation 



^sup V = c'^^'P 



ft V = c{x) 



flinf ^ — (^inf 



Figure 1: The propagation medium 

the Differential Calculus, sucli as: dt denotes tlie partial derivative in the time 
variable, if di stands for the partial derivative with respect to i = x, y, z, divv = 
dxVx + dyVy + dzVz rcprescuts the divergence of a vector field v with components 
VxiVy, Vz and the operator V is the gradient defined as Vi^ —* {dx'P, dycp, dz^p)- 
Model ([T|) can be derived from two laws of the mechanic of continuous media 
(see for instance [H]) which involve the velocity in the fluid u = u{t,x) =* 
{ux,Uy,Uz). The first concerns the conservation of the displacement quantity 
and reads as: for i = x,y, z, 

p r. = diva-i + }i (2) 

where / =* {fx,fy, fz) designates the density of the exterior forces and p rep- 
resents the density of the medium supposed to be constant herein. The stress 
tensor cr describes the behavior laws of the fluid whose lines are denoted by cr^ 
with i = x,y, z. We use the notation: 

,. daix , daiy daiz ,„, 

diver, = — l-^-^ + -s— ■ (3) 

ox ay oz 

In the particular case of a perfect fluid, the tensor cr can be written as : 

a- = -ph 
where I3 denotes the identity matrix of order 3. The law ^ modifies then into 

P^^-Vp + f. (4) 
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In the framework of seismic prospection, the exterior forces are often neghgi- 
ble. The source is an impulsion which is taken into account in the second law 
describing the conservation of the mass: 



1 dp 
v'^p dt 



^^ , divu = q, (5) 



Ht-t.y 



where q = q{t, x) represents an injection rate per mass unity. By applying the 
divergence to Eq.(|3]) and by deriving Eq.® with respect to the time, one can 
eliminate the velocity u and the acoustic pressure p is solution to the second- 
order wave equation with right-hand-side (rhs) given by S" = dtq. Thus, we 
have: 

q{t,x) = 5{x-Xs)^^ [e " 

The wave equation ([1]) can then be rewritten as a first-order hyperbolic system 

of the form: 

pdt (m) + Vp = dans ]0, T[x]R^, 

:^dtp + divu^q da.ns]0,T[x]R^, (6) 

p{0, a;) = et m(0, a;) = dans JR^. 

For the numerical simulations, the direction of the depth is selected as the one 
governing the sense of propagation. This approach is classical in case of homoge- 
neous media where the physical parameters do not vary. By using the formalism 
of pseudo-differential operators, it can be generalized to inhomogeneous media. 
In this work, we are interested in system ^ to which we associate a reduced 
system which is the subject of the following section. 

3 Reduced system 

The initial system is given by ^. We recall that the initial data are supposed 
to be null. We propose to describe the propagation of waves along the direction 
Oz. This idea was formerly exploited in [6| next in [3 [I] where the problem 
was to compute solutions propagating into plane and homogeneous layers. In 
such a way, the problem turned into the solution of a system with constant 
coefficients and the solution was written as a linear combination of eigenvectors 
of the problem. This approach was successful for the study of the propagation 
of plane waves into ID homogeneous layers [B] and the theory developed next in 
[7] provides an extension to the 2D case. Nevertheless, in case of more complex 
media, the constitutive parameters p and c of the domain f2 vary in all the di- 
rections. This is why we need the formalism of pseudo-differential operators to 
generalize this method of modelling, just as it was formerly suggested by M.V. 
de Hoop in [S]. 

In the following, we limit our attention to the solution in fl. We eliminate the 
time variable by applying a Laplace transform to the system and the related 
variable to t is denoted by uj. System ^ then transforms into a stationary 
system, which reads in Q as: 

'V ±p + ito (pu) =0, 

dzp + iu}{pul) = 0, 

^p + div± {u)' + dzU^^ q, 
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where V^ —* {dx,dy), u =* iux,Uy) and div±u = dxU^ + dyUy. In the above 
equations, ^ represents the Laplace transform of (p. The time variable being 
suppressed, the transverse unknowns u are eliminated by plugging the first 
equation into the third. Then we get a system with a pair of unknowns of the 
form: 

{D,+L)U^F (7) 

with U =* (p, v'z), Dz = ^2dz and I2 represents the 2x2 identity. The operator 
L is defined by: 




L 



and the source F is given by F =* (0, q) . 

By factoring by iw, L reads as L = zwL" where 



UJ 



is defined by: 







div I 





According to the theory developed by Hormander [T^] , Operator L^ is a pseudo- 
differential operator of OPS" depending on the parameter w, and if £" denotes 
its symbol, we have the following representation: for any test-function <p, 

iV = ^ \ \ O {x\ ^ V {s') e-<^' -'')■''' ds'dk' 

where s' =* {sx,Sy) e M and k' —^ {kx,ky) is the dual variable of x' such as 
the symbol of V^ is given by ik . 
Symbol of L* is defined by: 



with I fe' I — k^ + ky. Each term of symbol £" is an homogeneous function of 
k' with degree 0. Hence it belongs to a class of matrices whose terms are in 
S'^. Then it admits an asymptotic development with respect to the parameter 

\k' 

- — - which is small at high frequency, that is when w is large. The operator 
O can then be seen as a semi-classical operator. Moreover the symbol of L = 
iojL' admits an asymptotic development according to the small parameter oj~^. 
Hence we are in the framework of the w-calculus developed by O. Lafitte in [H] 
and L can be considered as a w-pseudo-differential. This property is interesting 
for the construction of high-order models and is the subject of a current work 

To solve (O turns into describing the propagation of the field U along the depth 
z. The formalism of pseudo-differential operators is well-adapted to extend the 
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classical way for solving differential systems with constant coefficients and based 
on the diagonalization of matrix L. In the case of pseudo-differential operators 
(or differential operators with variable coefficients ), Taylor [l^ has developed a 
method of factoring strictly hyperbolic systems. This method allows to replace 
the initial system by a system of two one-way equations. One will describe the 
down-going propagation and the other will model the up-going displacement. 
Both of the equations are coupled by a off-diagonal matrix which are assimilated 
to reflection terms. The idea of Taylor has been developed in [H [9], later in 
[15] for acoustic waves, in the simplest case where p is constant. Notice also 
that it has been applied in [2j for the construction of radiation conditions for 
electromagnetism. 



4 First-order Formulation 

The first-order formulation reads as a transport equation whose derivation is 
based on the factoring of the symbol £" . The eigenvalues of £" are given by the 
symbols 70 and —70 with 

In order to define the square-root involved in ([S]), the plane M^ is divided into 
three regions. The first is defined by 

■H= (fc'e^Mfc'l^ < 



c^ [x] 



As soon as fc' e H, matrix Cq admits two eigenvalues which are single and real. 
This region corresponds to the one in which System ^ is strictly hyperbolic 
and then Ti is called the hyperbolic region. It is exactly the region in which 
propagation occurs and the sign of each eigenvalue indicates the sense of prop- 
agation. The eigenvalue 70 is associated to the downgoing part of the wave 
i.e. the part propagating in the direction of increasing z. In the same way the 
eigenvalue —70 is related to the upgoing part of the wave. The second region is 
defined by 

£ = \k' eIR^,\k'? > 



(? {x) 



If fc' G 5, the eigenvalues of C\ are purely imaginary. System © is then elliptic 
in this region which is called the elliptic region. It defines a subset of frequen- 
cies which are not linked to the propagation. At last the third region is defined 

r, ,,2 iJ- \ _ 

\ys Q = <\k\ = -T—. — r > in which ([6|) degenerates. Indeed the eigenvalue is 
I ' c2(a;)J ^ 

double and the problem is neither hyperbolic nor elliptic. Region Q is called the 

glancing region and is related to grazing rays. 

Then we have the following result: 
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Proposition 4.1. Let rn G Z and U be a solution to the reduced system. There 
exists at least an operator Pm G OPS^^ with inverse Q-m G OPS^"^ such that, 
if V = Q-mU, V is solution to: 

{D,+iu;Ao)V^RoV + Q-^,nF (9) 

where Aq G OPS^ is a unique diagonal operator and Rq G OPS^ depends on 
p 

Proof. Let A^o be the matrix defined as 

^70 



-70 

We aim at constructing a matrix V such that tVJq = V^^C'^V. If w —^ (w+, w^) 
denotes one of the eigenvectors, associated with the eigenvalue ±70, its coordi- 
nates satisfy the equation: 

pw^ — ±7ow^ 

which admits an infinite number of solutions. We propose to solve this equation 
in such a way that the coordinates of the eigenvectors are in the same class of 
symbols. Hence if S"^,m G Z denotes this class [T^, we define the conjugated 
matrix Vm = ['U'"'',w~] where u""*" and w~ have been chosen in S"™. Then 
Vm defines an operator of OPS™ denoted by Pm- Matrix Vm is invertible 
and its inverse V^^ is the principal symbol of the inverse of Pm denoted by 
Q-m G OPS"™. As far as the symbols are concerned, we have the relation 

Vm^C^Vm = Mo- 

Then we define the operator Aq via its symbol by the relation: 

a(Ao)-A^o-'Tp(Q-m^*Pm). 

According to the pseudo-differential theory [iTj, we deduce that there exists a 
regularizing operator R i G OPS~^ such that 

Q_„,L«P,„ = Ao + ii-i. (10) 

Operator i?_i is uniquely defined by its symbol whose expression is known 
thanks to the composition rule of pseudo-differential operators [I^ . 
Next let us consider [/ as a solution to jT]). We set V = Q-mU. Then V 
satisfies: 

{D,+L)PmV^F 

and by composing on the left by Q-m, we get : 

Q^m{D, + L)PmV = Q-raF. 

Moreover since 

^-' z^ ra — -* m-'^z ^ ^z^ m 

where dzPm, is the pseudo-differential operator with symbol dzVm, the system 
modifies into 

{Dz + Q-mLPm)V + Q-mdzPmV = Q^mF. 
RR n° 0123456789 
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But by definition of L, we have : 

which imphes, by taking (fTO| into account, 

(D^ + iLoka + iojR-i + Q-md,P„r)V = Q-mF. 

Proposition 2.1 is then proved by setting : 

i?0 = —iuj{R^i) — Q^rndzPrn- 



Remark 1. Let m be fixed. Then, Matrix Vm is not single. That means that 
even if m is fixed, we can construct an infinite number of models which differ 
all by the coupling operator Rq. Nevertheless, it can be chosen in such a way 
the computational cost is a bit cut down, as seen further. 

Hence the model is defined from a transport equation which involves two 
one-way equations of order +1. In the rhs, Rq takes the coupling terms into 
account and they are modelled by the off-diagonal terms of the rhs. 
Our aim is to account for the lateral variations of the velocity. Hence the density 
mass can be constant or variable in all the directions. In this work, we limit our 
attention to the case p is constant. 

The coupling operator can be split into the sum of two operators, the first 
being defined by the diagonal part of Rq and the second by the off-diagonal 
part. Introduce 

R^ =: diag{Ro) and Hg'^ = Rq Rq- 

By considering this decomposition, one can construct a new model in which the 
transport equation is diagonal and involves the sum of two operators in OPS^ 
and OPS° respectively. Then we get: 

Corollary 4.2. Let m E Z. There exists at least an operator P,,,, G OPS™ 
with inverse Q-m G OPS^™ such that if U denotes a solution to the reduced 
system, V — Q-mU is solution to: 

{Dz + iLoKQ-R'^)V = RtV + Q-.„,F (11) 



The model in 14.21 consists of two one-way equations which are coupled by 
the operator R^'^. According to its definition, Rq is not easy to define exactly 
because it is given by an asymptotic expansion expressing its symbol. It is 
obvious that we cannot use the exact symbol of Rq and that we must be content 
with approximating the system by a truncated one. The truncation order must 
then be defined. Since the initial reduced system is defined from an operator L 
in OPS^ with a rhs in OPS*^, we propose to keep this property in the one-way 
model. Hence we address the question of putting the principal part iio,p in 
place of Rq. Operator i?o,p is defined as the one whose symbol is <jp{Rq). 
Then we get the following result: 
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Proposition 4.3. Under the same assumptions than in Proposition \4.1\ an 
approximate one-way model is given by: if V = Q-„JJ. V is solution to 

{D, + iwAo - R^p) V = W^%V + Q-^F, (12) 

where the symbol of Ro,p is given by: 

a{Ro^p) = -iu {v-J (v^./:») (v^.T'™) + (V;,.^;') v^^ (/:«n„) } - n;'a,K 

Proof. We aim at computing tlie principal symbol of Rq. We use the definition: 

Ro = -iuj (fl_i - (5_„i) d:,Pra 

with R I = Q -mL'* Pm ~ Ao, which implies that: 

a{Ra^p) = -iujap (-R-1+) - <Jp {Q-mdzPm) ■ (13) 

We introduce the following notation. Let A and B be two matrices whose terms 
are symbols. Let a —* {ax,ay) be a multi-index in IN"^. We set ^V^'AViB 
the product defined by: 



k X L — ^ k X 



l/3| = |a| 

where /3 =* {P^, (3y)GlN^, af, = 5f;af^'' and 9f, = 9^af . 

By applying the composition rule of pseudo-differential operators |17| . we get: 

ap {Q-rnd,Pr,r) = ^7/5.^.- (W) 

Moreover we have 

where AA-2 is the symbol of an operator in OPS^^. Then, since 

p-/£«n„ = ^(Ao), 
we can deduce that 

ap(i?_i) = ap(p„-i/:«n„-Ao) = -iv-,' {v ^, C^) {V^,P,^)-iV^,V,-'V^> {C^Vm) 

(15) 

By plugging (fTSJ) and (fT4|) into (fT3|) . we complete the proof of Proposition 

1431 ■ 



5 Setting of the numerical method 

The numerical method is now based on the solution to the approximate one-way 
model: 

{D, -f luAo - R^p) V ^ {Rip + R^%) V + S (16) 

where Aq is the diagonal operator in OPS° whose symbol is the diagonal matrix 

-»^(? -V 
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Symbol 70 = y ^ — -^ and its opposite —70 are the eigenvalues of £'. The 
auxiliary unknown V whose components are the down-going field V+ and the 
up-going one V^ is linked to U by the relation U = PqV where Pq is the 
operator in OPS° whose symbol is the matrix Vq. We choose Vq such as 

Kp = ~\^o^d,TQh and R'^'^p - ]^T^,^d,ToJ2 

where Fq G OPS" is defined by symbol 70. Then, Vo is given by "Pq = 

1 . Since R^ p does not play a part in coupling the down-going 

and up-going fields, it can be seen as a transmission term. Introducing a pa- 
rameter e G {Oi 1}: we write 

(£), + iLoka - eRIp) y = ((1 - e)Rlp + Rl%) V + S. (17) 

Then if e = 0, Rq p acts in the right-hand-side of the model (like in the model 
proposed by De Hoop), and if e = 1, Up p is included in the left-hand-side (like 
in the one-way equations proposed by Zhang et at). Only Rq'^p accounts for 
the coupling and it can be seen as the reflection operator. In the following, we 
use the notations: 

jyd T' T ryad T? J 

-"-O.P ~ -'0-'2 , -Kg p — UqJ2 

with the letters T to indicate the Transmission and R corresponds to the Re- 
flection. 
We introduce also the operator: 

AP = iujAo - eTo/2- 

To compute the solution to (fTC|) requires to invert D^ + M^ and the inverse 
operator is denoted by G*^, the so-called propagator. It governs the propagation 
along the depth, in the two senses. Like M*^, G^ is diagonal and 

Assume that C^ is given. Then p6p can be transformed into: 

V^G' ((1 - e)nh + R0J2) V + G'S 

and if we introduce K'^ ^ C^ ((1 — e)To/2 + R0J2), we have 

{I2 - K')V = G'S. (19) 

Hence the solution is computed once {I2 ~ K'^) has been inverted. By construc- 
tion, K'^ £ OPS~^ because K'^ arises from the composition of G*^ G OPS^^ and 
of ((1 - e)Tf)l2 + R0J2) e OPS°. Hence we have, according to [H]: 

il2-KT'=J2(^y- (20) 

j>0 
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This representation corresponds to a Neumann series for the inverse of I2 — K'^ 

under the assumption \\K'^\\ < 1. 

Then we use (|20|) to expand Y in the foUowing form: 



{ j>o (21) 

i v" = G''S et y(j') = 7^^y(j'-i),j > 1 

This representation of V is called the Bremmer series and refers to as pre- 
cursory Bremmer's works [5] in the simple case of a ID model. 
Each term V*^^^ has two components denoted by V_^ and V_ where + corre- 
sponds to the down-going part propagating in the sense of increasing z and — 
is associated to the up-going part. 

By definition of the source F =* (0,g), the auxiliary source S has a priori two 
components 5*+ et 5_ and the first term in the series is given by: 

K[°^ = G;5+ et Ki°) = CLS^. 

Then the second term V*^^^ is obtained by the formulas: 

where Kf^ , 1 < l,m < 2 represent the terms of K'^ . 

However the numerical solution is computed into a region limited by z = and 
z = Zraax and which is surrounded by two homogeneous regions. By definition 
V_ represents the propagation of the component S^ of the source which is 
located at the surface z — Q. Hence the support of V_ is embedded in ^sup- 
In the numerical tests we present, we do not evaluate V_ . The algorithm is 

evaluated with V_ = which does not produce any error in the numerical 
results because the receivers are located at the same depth than the source. 
Nethertheless the omputaion of V_ is possible and the case of a source located 
into f2 can be considered also. If K'^ is represented from its principal symbol: 

ap{K') = ap{G') {a{To){l - €)h + fT(i?o) J2) , 

we have ap{K'^) = in Qgup- Then since K""- is a pseudo-differential operator 

depending continuously on the parameter z, iff2^- 0: 1 < ' < 2. We then 

deduce that in that case, it is not useful to compute V_ . Then V'^^^ reads in 
the simplified form as: 

vi'^^Kl.vi'^ and vl'^ ^ K^,,vf 

while the next terms are linear combinations of the down-going and up-going 
fields: 

with j > 2. If e = 1, each term is simpler because Kf^ = K22 = 0, which implies 
that V_|: — and next, using the chain rule, we get: 

Vi''+'^ = et V}"'^ - , J > 0. 
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Hence each term is defined from either a down-going field or a up-going one 
which shows off the uncouphng of the computations. This property makes 
think out models based on paraxial equations and a comparison between one- 
way models and paraxial systems is a current work [4J . 

To summarize, the solution of P^ involves the following steps: 

1. decomposition of the source: S = P^^^F 

2. computation of the Brcmmer terms from the propagator G"^, the reflection 
operator Rq and the transmission term Tq. 

3. recomposition of the solution: U = PqV. 

The center of the algorithm is given by the item [21 In order to illustrate its prin- 
ciple, we consider the case of a velocity model consisting of three homogeneous 
layers. 

6 Bremmer series for a 2D stratified medium 

The numerical solution is based on the expansion of the fields as a Brcmmer 
series. We assume that ft is defined as in Fig. [2j The first layer is homoge- 
neous with thickness zi — ^ where J is a small parameter and the corresponding 
constant velocity is denoted by ci. Next on a thin layer with thickness (5, the 
velocity continuously varies from ci to the constant value C2 which is the prop- 
agation velocity in the layer with thickness {z2 — zi) — S. Then on a thin layer 
with thickness S, the velocity continuously varies from C2 to the constant value 

C3- 

The interfaces between each medium are fiat and they are linked to varia- 
tions of the velocity which arc so important as S is small. 

As long as the wave propagates into a homogeneous layer, no refiection phe- 
nomenon occurs. This is well-reproduced by the symbol Rq since it is equal 
to as soon as the z-derivative of the velocity vanishes, as in the case of a 
homogeneous medium. Then we have: 

r A ^ r A A"i r A ^ 

a{Ro) = a{To) = if z e }^z < z^ - -j\J^z, + -, Z2 - -^\J ^z > z^ + -j 

(22) 

6.1 Series terms before any discretization 

We propose to write the terms of the Bremmer series with item 0, 1 and 2. The 
next terms can be straightforwardly deduced from the term with item 2. Let us 
assume that if ip =* {'P+, (/?-) is a test-function depending on x' and z, we have: 

G;(^+= / G%{h, z)ip+{h)dh (23) 

Jo 

and 

Gtip- ^ [ G'_{h, z)ip^{h)dh (24) 
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Figure 2: Velocity profile for ci < C2 < C3 



where any G'j.{h, z) is a pseudo-differential operator depending continuously on 
h and z. It has been proven in ,16J that this assumption is available. Then 
the first term in the Bremmer series is obtained from the following arguments. 
First, just as was formerly noted, it is not useful to compute V_ . This is why 
we set: 



Next according to its definition, Vi reads as: 



VO < z < z. 



.(0) 



max 5 * -\- 



vr{z) 



G%{h,z)S+ih)dh 



As far as the second term of the series is concerned, we have: 

VO < z < z^a. , vi'Hz) = (i^l^y|"^) (z) , Vi'\z) = [k^.vI"^) (z) 



where: 

Kl^ = (1 - e)G%To and K^^ = CLRq. 

Firstly consider the down-going term V_^ . Then, Property 
(ToVI^A (z) = for all z e [0, zi - 6/2] . 
and we deduce from (1^ that: 



implies that: 
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V^'^ (z) = for < ^ < zi - 5/2 
which gives rise to 



Vl^'iz) = (1 - 6) / GXih,z){ToVl''>){h)dh 

Jzi-5/2 

if z G [zi — (5/2 < z < zi + (5/2]. Next using (221) again, we get: 
(ToVI"A (z) = for all z e [zi + (5/2, za - (5/2] 
and combining this relation with ()23p yieMs: 



pZi+5/2 

Vi'\z) = (1 - e) / G;(/i,z)(roF|"^)(/.)d/i 

■/zi-(S/2 

for any z into the homogeneous layer [zi + ^/2, Z2 — (5/2]. 
Now, if z belongs to [z2 — 5/2 < z < Z2 + 5/2], we have: 

pzi+S/2 

Vi'\z) = (1-e) Gl{h,z){nvi°^){h)dh 

Jzi-S/2 

+ (1-e) / G%{h,z){nvl°'>){h)dh 



IZ2-S/2 

and in the homogeneous layer {z > Z2 + 5/2} 

rzi+5/2 

Vi'\z) = (1-e)/ G%{h,z){ToVi'>^){h)dh 

Jzi~S/2 
nZ2+S/2 

+ (1-e)/ GX{h,z){ToVf){h)dh 

J Z2-S/2 



If we choose t — Q,V/^ corrects V"^ in which only the downward propaga- 
tion of the source is taken into account. Hence v!^ plays the role of a corrector 
of the propagation by introducing the transmission effects. 



Next when e—\,V\_ vanishes. In fact the transmission is included into the 
propagator G\ and V]^ directly accounts for the transmission effects at each 
interface. 
Now consider the up-going term. According to ([M|) . 

V[^\z) = / G'_{h,z){R^V^^){h)dh ^ „ /■ ""^ G'_{h,z){RoV^f){h)dh. 



Moreover according to (P^ , we have: 

V[^\z) = for z^ax >z>Z2 + 5/2. 
By applying the same reasoning than previously, we get: 

V^i\z) = - Gt{h,z){RoVl°^){h)dh forz e [z2 - 5/2, Z2 + 5/2], 
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V^^\z)^^ / G'L{h,z)RoVl°\h)dh ior z e [zi + S/2,Z2 - 6/2], 

JZ2-S/2 
J-Z2+S/2 

V^^\z) = - / Gt{h,z){RQV'^^){h)dh 

Jz2-S/2 
zi+S/2 

G'L{h,z){RoVl°^){h)dh for z e [zi - 6/2,zi + S/2], 

J z 

and 

PZ2+S/2 

Vl^\z) = - / Gt{h,z){RoVl"^)ih)dh 

Jz2-S/2 
zi+S/2 

G'Lih, z){RnVl"'){h)dh for z e [0, zi - (5/2]. 

lzi-5/2 

In the case where e = 0, G" represents the propagation without accounting 
for the discontinuity of the medium. Hence V_ results from the propagation 
of the reflected part oi V_\_ . 
To construct V'^^ we follow the same scheme than for V^^'^: 



Vl'^ ^ Kl.vi'^ + KI,V1'\ 



(2) 



The component V}_ reads as follows: 



VO < z < Zrna. , Vl^'iz) - / G%ih,z)iil - e){nVl''){h) + {R„Vl'')ih))dh. 

Jo 

In the first layer, we have seen that cr(To) = cr(_Ro) = 0. Thus we can deduce 
that 

V0< z<zi -(5/2, Vl^^z) ^0. 

Next developing the same ideas than for the construction oi V^ , we get the 
following relations: 

for zi-S/2 <z < zi+ (5/2 

Vi'\z)^ r GXiKz)({l-e)iToVl'\h) + iRoV}_'^){h))dh, 



/zi-S/2 

for zi+S/2 <z < Z2- (5/2 
rzi+S/2 

I zi-S/2 

for Z2 - <5/2 < z < Z2 + (5/2 



Vi^\z) = r G^iKz) (^{1 - e)iToVi'^)ih) + (i?o^i'^)(M) dh, 



Vi'\z) = r G%{h,z)(^{l-e){ToVl'^){h) + {RoV^'^)ih))dh 



G%ih,z) (a - e)iToVi}^)ih) + iRoV^'^){h)) dh, 

Z2-S/2 ^ ' 
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and for z > Z2 + S/2 

fZi+S/2 

Vl^\z) = / GX{h,z){RoVl^^){h)dh 

Jzi-5/2 
/•z+S/2 

+ / G%{h, z) ( (1 - e){nvi'^){h) + (i?oVd'^)(M ) dh. 

Jz2-S/2 ^ ' 

(2) 

Formulas describing V\_ into the different layers are more intricate than the 
ones given V\_ because of V_ ^0. Now the down-going term is a linear com- 
bination of down-going and up-going terms. 

In the same way, V_ ' satisfies the relation: 

VO < z < z^a. , Vl^\z) ^ _ /■ ""^ G'_{h,z){{l-e)(ToVi'^){h)+(RoVl'^){h))dh 



which can be split into each layer as: 

if Zmax > Z > Z2 + 6/2 



Vi^\z)^0, 



if Z2 + 5/2>z>Z2- s/2 

rZ2+S/2 



Vl'\z) = - / GtiKz) ((1 - e)iToVl'^)ih) + (i?oK['')(M) dh, 



rZ2+d/2 

Vi'\z) - - / GUh,z) [^{1 - e)(ToVi'^)(h) + (RoVi'^){h)) dh, 



if Z2-S/2>z>zi+ 5/2 

rZ2+S/2 
lz2-S/2 

if zi+5/2> z>zi- 5/2 

rz2+S/2 



V}_^\z) = - Gtih,z)Ul-e)iToVP)ih) + iRoVi}^){h))dh 

Jz2-S/2 ^ ^ 

pzi+S/2 

-J GUh,z) [^{l - e){ToV['^)(h) + {RoVi'^){h)) dh 

and if zi + (5/2 > z > 

V}_^\z) = - r'^"Gtih,z)(^{l-e)iToVl'^)ih) + iRoV}'^){h))dh 



Z2+S/2 

Z2-S/2 
l+<5/2 

GUh,z) ((1 - e){ToVi'^)(h) + (i?o^|'^)(M) dh. 

zi-S/2 ^ ' 

In the case where e = 1, V_: vanishes since V\ vanishes too. 
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The next terms V^ and V_ with j > 3 read exactly in the same way than 
in the case j = 2 but by replacing the subscript (2) by (j) and the subscript (1) 
by (j — 1) into the formulas. 

In the case where e — 1, the method requires two times less operations than 
in the case where e = 0. Hence one could already suppose that it will be more 
judicious to choose e = 1. This is why we propose to investigate this question 
in the following section devoted to numerical experiments. 

6.2 Numerical Bremmer terms 

Let Az be the depth step. In practise, Az — 6 where 6 is the thickness of 
the layer in which the velocity varies continuously between two constant values. 
Each of the integrals involved into the Bremmer terms is approximated by the 
rectangle method and to choose Az = S implies that the integrals in the variable 
z is replaced by the multiplication by Az. 

For instance, when the receivers are located at the surface z = 0, they are 
supposed to record V-(0) which is obtained step by step by computing the 
terms of the Bremmer series and we propose to truncate the series up to the 
third term. The approximate field is then given by: 



with 



(0), 



Vl^\0) = -Az(Gt{z2,0)RoVi°\z2) + Gt{zi,0)RoVi°\z,)y 
Vl^\0) = -Az(Gt{z2,0)RoVi'\z2) + {l-e)Gt{z,,0)ToV['\zi)y 
v[^\q) = -Az(Gl(z2,0)i?oV:['^(^2) + (l-e)Gl(zi,0)ToT/i'^(zi)). 

The terms are obtained by following always the same scheme. We focus our dis- 
cussion on V_ computed at a given depth z. Nevertheless the scheme depicted 
at Fig. [3] gives a general survey of the contribution of each term of the series. 
Assuming the propagator Gj-(., .) satisfies the Chasles relation: 

G'^{z,h)G'^{h,z*) = G%{z,z*) = G%{h,z*)G'±{z,h), (25) 

the approximate term can be written as a function of the source: 

y_,app(0) = -Az(Gl(zi,0)i?oG;(0,zi)5+) 

^ ^ -f 

Reflection at the first interface 

- Az (Gl(zi, 0)(/d - Az(l - e)To)GUz2, zi)RoG%{zi, Z2){Id + Az(l - e))ToG;(0, zi)S+) 

Reflection at the second interface 

+ (Az)3Gl(z2,0)i?oG;(zi,z2)i?oGl(z2,zi)i?oG;(0,z2)5+ 
^ ^ ^ 

First multiple 
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Figure 3: The Bremmer terms 
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In Geophysics, the two first terms are the most important. But it is also of great 
outstanding to compute the multiples also because in the case where the velocity 
varies strongly, they can have a large amplitude and a simultaneous calculus of 
the two first events with the multiples can give rise to polluted results which are 
very difficult to interpret correctly. Hence it is very interesting to dispose of a 
numerical method which allows to compute the multiples and to uncouple them 
from the first events. Then the multiples can be removed from the primary fields 
which facilitates the interpretation of the results. To compute the multiples, it 
is necessary to combine the reflection operator at least three times. Precisely 
the TO-th multiple is the result of 2m + 1 reflections. 

In order to explain the formula giving V^^app{^)^ we restrict ourselves to write 
V_ {z) at the depth z. The other terms are obtained in the same way. Before 
discretization, V_ (z) is given by: 

V[^\z)= r Gl{h,z){RoVl"^){h)dh. 

According to the Chasles ([25ll relation, we have: 

Vz*e[z,z™aJ, Gl{h,z)^Gl{z*,z)Gl{h,z*) 
which implies 

V^^\z) = G'_iz*,z) r G'_{h,z*){RaV^^){h)dh+ C G'_{h,z){RoV'f^){h)dh 

Thus we get the relation: for any z* G [z, z^ax], 

V^^\z) = G'_{z* , z)v!:^\z*) + I G'_{h,z){RoV^^){h)dh 

J z* 

In practise, z* — z + Az. Hence if V_ {z) denotes the discrete value of Vj; , 
we have by approximating the second integral by a rectangle method, 

V^^\z) = Gl{z + Az, z)v'l^\z + Az) - AzGi(z + Az,z){RoVl°^){z + Az). 



The case of a stratified medium is very useful to show clearly why the cases 
e = and e = 1 are different. Indeed, since the symbols do not depend on the 
variable a;', we can write the propagator in the wave number domain and the 
symbol to of Tq is non zero only at each interface. Hence if we consider a layer 
{zi < z < Z2}, the propagators G^. are represented in the Fourier variables as 

[simiii]: 



and 



„c _ g-i(z2-zi)tJ70geAzto(2l) 



„« _ g-i(z2-zi)ct;70g-£Azto(zi) 



By injecting these relations into the definition of V^-,app(0), we get for the pri- 
mary refiections: 

• when e = 
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• and when e = 1 

^ g-iziW70g-Azto(2i)g-i(z2-zi)'J7o/'_A2„ C~ '\~jg-J(z2-zi)w70g-Azto(zi)g-iziW7o3 



We can then observe that the term (/d=F (^^io(-Zi))) acting in the case e = is 
replaced by the exponential eTAzto(zi) jj^ ^j^g ^g^^ e = 1. These two terms are 
of the same order when l^ztQ{zi) is small, which is the case when the medium 
is smooth, i.e. when the vertical variations of the velocity are small. 

7 Description of the software 

We use the language Fortran 90. The code has been developed from the C 
one written by Jerome Le Rousseau. It runs in the 2D and 3D cases. Its 
architecture has been modified as compared to the first version. Each routine is 
now independent of the other ones and can be replaced by anyone else. This is 
an interesting property for the analysis of the numerical method which can now 
be coupled with other methods like the ones involving finite difference schemes 
for instance. The computation of the Bremmer series has been optimized in the 
case where e = [16]. Hence, the computational cost is almost the same for 
e = and e = 1. 

Data 

The source F is obtained from a FFT (Fast Fourier Transform) of the source 
q. One gets a table as a function of the frequency uj. The source g is a Ricker 
function which implies that its FFT decreases quickly to zero. This is why we 
define a window of calculus into the interval [0, Umax]- Next one sets the number 
of terms in the Bremmer series by pointing at the number N of multiples which 
must be modelled. Then the code will compute V_ up to j — 2N + 1. 

Decomposition of F 

The cases e = and e = 1 are distinguished. In both cases, oj varies from to 
uj„iax and 

(i) e = 1 We know that when e = 1, T/|^^+^) = and ^1^^' = 0. 

The integer m varies from to TV and one computes V_l_ et V_; *" at 
each depth llS.z where I varies from to Imax with Zmax = Imax^z. This 
can be summarized by: 

Initialization 

for / varying from 1 to Imax 
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compute V^\ll\z) = G\ {S+{{1 - l)Az)) 

store RqV'^^ 
stop the variations of I 
for / varying from Imax — 1 to 

compute V[^\l/:^z) = Gl (RoVI"'' {{I + l)Az) 

store RoV[^^ 
stop the variations of I 
for m varying from 1 to N, 



for I varying from 1 to Imax 








compute y|^"'(;Az) = 


ov(. 


RoKi'"" 


-'Hii 


store RoVl_ ™ 








stop the variations of I 








for I varying from l^ax — 1 to 






compute yi""+'^(/Az) 


= Gl 


[Rovr 


-\{i 


store i?oVd'"+'' 









l)Az) 



l)Az) 



stop the variations of I 

The results are stored in a table RESULT ■- RESULT{x', m, uj). 
stop the variations of m 

A priori it is the most complicated case because the downward and 



upward parts of each term are coupled. However the algorithm proposed 
initially by J. Le Rousseau can be optimized. The optimization is based on 
a summation process allowing to compute both the 2j-th and the 2j + 1-th 
terms in the same time. This reduces considerably the computational time 
when e = and makes it competitive with the case e = 1. 

Computation of the Bremmer terms 

Let j be the summation index into the series. We use the temporary table 
TEMP := TEMP{x') and the results are stored into the tables TAB+ 
and TAB^ with length Imax- 
for m varying from to A^ 

Initialization of the downgoing part 
TAEMP := 
TAB_(0) :=0 

for ? == 1 to Imax, 

if TO = 0, TAB+{1-1) :=0 

TEMP := G°_^{{l-l)Az,lAz){TEMP-TAB^{l-l)+TAB+{l- 

1)) 

TAB^{1) := -AzRqTAMP 

End of variations of I 
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Initialization of the upgoing part 

TEMP := 

TAB+{l„^ax) :=0 
for I = l„iax to 1, 

TAB+{1) ■- AzRqTAMP 

TEMP := G°_{lAz, {I - l)Az){TEMP ~ TAB+(l) + TAB^{1)) 

End of variation of / 

The final values of TEMP give the values of the m-th multiple which 
are stored in the table RESULT := RESULT{x',m,uj). 

End of variations of m 

Composition 

For each value of to, we have computed the field recorded at z = and this has 
been done for each value of the frequency from to Umax- Then applying Pq 
to this field, we get the unknown U. Next by applying an inverse EFT to its 
first component, we get the value of the acoustic pressure as a function of time. 
By depicting the acoustic pressure at the surface z = we get a seismic section 
which represents the time value of the pressure as a function of the transverse 
variable x' . In all the pictures we state the variable x' is the abscissa and the 
ordinate gives the time oriented to the bottom. In this way we can represent 
the kinematic of the phenomenon. To depict the dynamic, we use a grey scale 
for the amplitude of the acoustic pressure at point {x' , t) and z = 0. 
It is also possible to represent snapshots (value of the wave field at (x',z,t)). 
Then we use an auxiliary table TEMP := TEMP{x',l,u!) in which we store 
the values of the field at {x', z, uj) and by applying an inverse EET, we get the 
field evaluated at {x',z,t). 

8 Setting of the velocity model 

We consider a simple velocity model (See EigH]) which consists of two layers. 
The topographic data of the model are the following: 



depth 


5 km 


thickness 


21.59 km 



The shot is set at the middle of the domain i.e. at Xg = (10800,0) and the 
receivers are located at the depth z = SOOOto. By this way, we only account for 
the transmission effects. The transverse dimension of the model has been chosen 
such that the cardinal of the discrete set of points Xj is matched for using FFTs. 
Each layer corresponds to a constant value of the velocity and their interface 
is located at z = 2500m. The mesh dimensions are collected in the following 
table: 
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[ FFT in time of the source .F | 



Loop on the frequences: operations representation by FFT 



Decomposition of the source ig 



Loop on the multiples number m 



Loop onl <l < Imax} 



Propagation from i — 1 to / of the downgoing field V} "^ 



Reflection: couplage with the upgoing field V_ 



(2m+l) 



I Update of the field with the reflection part computed at the step (m — 1) 



Transmission: couplage with the downgoing field V[ "^ 



End of the loop on I 



End of the downgoing part: champ= 



Loop on Imax > ^ > 0] 



Reflection: couplage with the downgoing field V 



(2m + 2) 



Transmission: couplage with the upgoing field V_ 



(2m + l) 



/< K 



Update of the field with the reflection part computed at the downgoing step 



Propagation from Z to Z — 1 of the upgoing field V_ "* 



End of the loop on I 



Record the data at the receivers 



End of loop on the multiples 



End of th loop on the frequencies 
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Figure 4: A two layers model. 
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width of the computational domain 


21.59 km 


transverse step Ax 


10 m 


depth step Az 


10 m 



In the foUowing we refer to Velocity Model 1 (VM 1) when the velocity is equal to 
1600 ms~^ in the first layer and 2400 ms^^ in the second one. VM 2 corresponds 
to the case where the velocity is equal to 1600 ms~^ in the first layer and 5000 
ms~^ in the second one. The last model VM 3 is defined by a velocity equal to 
2400 ms^^ in the first layer and 1600 ms^^ in the second one. VMs 1 and 3 are 
examples of velocities with a quite low contrast while VM 2 gives an example 
of high velocity contrast. 

On FigEl we have collected four scismograms. The first one depicts the 
reference solution obtained from the second-order full wave equation. On the 
right of this picture, we have set the solution of the one-way equations system 
as it was computed initially in [TS] . One can observe that the kinematics is well- 
reproduced. Nevertheless the extrema of the grey scale shows that the amplitude 
of the one-way solution is erroneous. This motivates the two following pictures 
underneath. The left seismogram depicts the best results: both the kinematics 
and the amplitude of the wave field are correct. The right seismogram is not so 
bad even if the amplitude is more erroneous than in the previous case. Hence 
this first collection of numerical tests shows that to include the transmission 
term into the model really improves the amplitude of the wave field. Moreover 
if we limit our analysis to consider seismograms, we can conclude that the 
transmission term can be numerically handled as a term of the right-hand-side 
of the one-way system or as a proper term of the one-way equations. But the 
seismograms give a global view of the results and they are not precise enough 
to estimate the accuracy of the amplitude. This is why we have performed a 
series of numerical tests and we have chosen to represent them from the value 
of Q which is defined by: 

maxt|pft,x,0)|^„„^„„^ 
maxt|p(i,a;,0)|„„^_^^j,' 

On Fig. m the lower curve depicts the value of Q for the one-way solution 
with e = 0. Its values are very close to 1 which confirms what was observed on 
the corresponding seismogram. The top curve represents the values of Q for the 
one-way solution without the transmission term. This picture strengthens the 
previous conclusion since the minimum value of the error is 20%. Hence it is 
essential to include the transmission into the model to get accurate amplitudes. 
Now, being convinced that the transmission must be included, the theory shows 
that it can be equivalently introduced whether into the one-way equations or 
into the rhs of the system. The first way allows to write a system of one-way 
equations which are quite close to the equations considered by Zhang and al. 
[18j as an approximation of the full wave equation. The second one is the most 
natural in the formalism of Bremmer series and was introduced by Le Rousseau 
and De Hoop [ISIE]. 

On Fig. [71 we have collected the values of Q for e = and e = 1. When 
e = 1 , the results are correct in a neighbourhood of the shot and then they spoil 
quickly. Hence this numerical test indicates that for VM 1, the best approach 
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Figure 6: Values of Q{x) in VM 1. Top curve: one-way solution without trans- 
mission. Lower curve: one-way solution with transmission (rhs). 
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Figure 7: Values of Q{x) for the one-way solution with transmission in VM 1. 
Top curve: e = 0. Lower curve: e = 1 
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Figure 8: Values of Q{x) for the one-way solution with transmission in VM 2. 
Top curve: e = 0. Lower curve: e = 1 

consists in including the transmission into the rhs. 

Thus the result for e — 1 may deteriorate more and more as fast as the velocity 
contrast is high. That motivates the next test on Fig. [5] where we can observe 
that actually the results spoil more and more: the minimum value of the error 
is now 8% versus 1% for VM 1. 

On Fig. [9l we have collected three curves obtained from the one-way solu- 
tion with e = 1. Each of them are distinguished from the value of the velocity 
contrast. The pattern shows that the lower the contrast is the larger the neig- 
bourhood of the shot in which the results are satisfactory is. This figure seems 
to indicate that the model with e = 1 is correct when the propagation medium 
is smooth. 

On Fig. [To] we have depicted the results for VM 3. It is just to show that 
the same conclusion holds even if the upper velocity is larger than the lower 
one. Let us mention that each curve shows some instabilities (depicted as local 
maxima) which are due to the periodicities created by the FFTs. 

9 conclusion 

In this paper we consider the numerical analysis of two one-way systems derived 
from the general modeling of M. De Hoop ^8^. 

Such a formulation is used to replace the full wave equation by a system of 
one-way wave equations whose computational burden is lower than the one 
associated to the finite difference solution of the wave equation. Moreover it 
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Figure 9: Values of Q{x) for the one-way solution with e 
velocity contrasts 
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permits to unpack and identify the multiples from the primary reflections. We 
have include a transmission term in the one-way model. Then numerical tests 
have been performed in the 2D case and they indicate that accounting for the 
transmission improves significantly the amplitude of the solution. The com- 
putational algorithm has been optimized in such a way that its complexity is 
now of the same order than the one of the one-way solver without transmission. 
Hence to add the correcting term does not penalize the computational method. 
Now the solution is based on a method of propagators whose numerical approx- 
imation accuracy is very sensitive to the position of this correcting term. To 
put it into the one-way system, which is the case e = 1, seems to deteriorate 
the results. But this approach should be interesting since it corresponds to the 
same idea than the one proposed by Zhang et al. in [TH] where time-arrivals 
are computed from the solution of a second-order wave equation obtained by 
factoring the full wave equation. 

In this paper, we have performed a numerical test which shows that the 
accuracy of the method with e = 1 is similar to the one of the method with 
e = when the velocity contrasts decrease. 

In the proposed methods, the transmission operator is approximated by a zero 
order pseudodifferential operator which is exact for stratified media. A higher 
order operator, that should account for media with lateral velocity variations, 
is the topic of a current research [3] . 
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Figure 10: Values of Q{x) for the one-way solution with transmission in VM 3. 
Top curve: e = 0. Lower curve: e = 1 
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